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Abstract 

A stochastic model for a chemical reaction network is embedded in a one-parameter 
family of models with species numbers and rate constants scaled by powers of the 
parameter. A systematic approach is developed for determining appropriate choices 
of the exponents that can be applied to large complex networks. When the scaling 
implies subnetworks have different time-scales, the subnetworks can be approximated 
separately providing insight into the behavior of the full network through the analysis 
of these lower dimensional approximations. 
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1 Introduction 



Chemical reaction networks in biological cells involve chemical species with vastly differing 
numbers of molecules and reactions with rate constants that also vary over several orders of 
magnitude. This wide variation in number and rate yield phenomena that evolve on very 
different time-scales. As in many other areas of application, these differing time-scales can 



be exploited to obtain simplifications of complex models. Papers by Rao and Arkin (2003) 



and 


Haseltine and Rawlings 


(2002 


) stimulated considerable interest in this approach and 


notable contributions by 


Cao, Gillespie, and Petzold 


(2005 


), 


Goutsias 


(2005 


), 


E, Liu, and 


Vanden-Eijnden 


(2007 


)• 


Vlastny, Haseltine, and Rawlings 


(2007 


), 


Crudu, Debussche, and 
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Radulescu (2009), and others. All of the cited work considers models of chemical reaction 



networks given by continuous time Markov chains where the state of the chain is an integer 
vector whose components give the numbers of molecules of each of the chemical species 
involved in the reaction. Most of the analysis carried out in this previous work is based 
on the chemical master equation (the Kolmogorov forward equation) determining the one- 
dimensional distributions of the process and is focused on simplifying simulation methods for 



the process. In contrast, the analysis in Ball, Kurtz, Popovic, and Rempala (2006), is based 



primarily on stochastic equations determining the process and focuses on the derivation of 
simplified models obtained as limits of rescaled versions of the original model. 

The present paper gives a systematic development of many of the ideas introduced in 



Ball et al. (2006). First, recognizing that the variation in time-scales is due both to variation 



in species number and to variation in rate constants, we normalize species numbers and rate 
constants by powers of a fixed constant N which we assume to be "large." 

Second, we replace N by a parameter N to obtain a one-parameter family of models and 
obtain our approximate models as rigorous limits as iV — > oo. It is natural to compare this 



approach to singular perturbation analysis of deterministic models (cf. Segel and Slemrod 



(1989)) and many of the same ideas and problems arise. This kind of analysis is implicit in 



some of the earlier work and is the basis for the work in 
Third, as in 



Ball, Kurtz, Popovic, and Rempala (2006), the different time-scales are iden- 



Ball et al. (2006) 



tified with powers Nq, and making a change of time variable (replacing t by tN 1 ) we get 
different limiting/approximate models involving different subsets of the chemical species. As 
observed in Cao, Gillespie, and Petzold (2005) and E, Liu, and Vanden-Eijnden (2007), the 



variables in the approximate models may correspond to linear combinations of species num- 
bers. We identify the time-scale of a species or a reaction with the exponent 7 for which the 
asymptotic behavior is nondegenerate, that is, the quantity has a nonconstant, well-behaved 
limit. The time-scale of a reaction is determined by the scaling of its rate constant and by 
the scaling of the species numbers of the species that determine the intensity/propensity 
function for the reaction. The time-scale of a species will depend both on the scaling of the 
intensity/propensity functions (the reaction time-scales) and on the scaling of the species 
number. It can happen that the scaling of a species number will need to be different for 
different time scales, and a species may appear in the limiting model for more than one of 
the time scales. 

Fourth, the limiting models may be stochastic, deterministic or "hybrid" involving stochas- 
tically driven differential equations, that is, piecewise deterministic Markov processes (see 



Davis (1993)). Haseltine and Rawlings (2002) obtain hybrid models and hybrid models have 



been used elsewhere in reaction network modeling (for example, Hensel, Rawlings, and Yin 



(2009), Zeiser, Franz, and Liebscher ( 2010[ ) ) and are a primary focus of Crudu, Debussche 
and Radulescu| ( |2009p. 

Finally, as in Ball et al. (2006), we carry out our analysis using stochastic equations of 
the form 



X(t) = X(0) 



k 



Y, 



\ k {X{s))ds)C k 



that determine the continuous time Markov chain model. Here the are independent unit 
Poisson processes and the (k are vectors in Z d . These equations are rescaled and the analysis 
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carried out exploiting the law of large numbers and martingale properties of the (For 



more information, see Kurtz (1977/78) and Ethier and Kurtz (1986), Chapter 11.) The 
other critical component of the analysis is averaging methods that date back at least to 
Khas'minskii| Ql966a|b"l ). (We follow [Kurti] ( p92| ) . See that paper for additional references.) 

If N is large but not large enough, the limiting model obtained by the procedure outlined 
above may have components that exhibit no fluctuation but corresponding to components 
in the original model that exhibit substantial fluctuation. This observation suggests the 
possibility of some kind of diffusion/Langevin approximation. Under what we will call the 
classical scaling (see Section |2j), diffusion/Langevin approximations can be determined simply 
by replacing the rescaled Poisson processes by their appropriate Brownian approximations. 
In systems with multiple time-scales that involve averaging fast components, fluctuations 
around averaged quantities may also contribute to the diffusion terms, and identifying an 
appropriate diffusion approximation becomes more delicate. These "higher order" correc- 



tions will be discussed in a later paper, Kang, Kurtz, and Popovic (2010). 



Section [2] introduces the general class of models to be considered and defines the scal- 
ing parameters used in our approach. For comparison purposes, we will also describe the 
"classical scaling" that leads to the deterministic law of mass action. Section [3] describes 
systematic approaches to the selection of the scaling parameters. Unfortunately, even with 
these methods there may be as much art as science in their selection, although perhaps we 
should claim that this is a "feature" (flexibility) rather than a "bug" (ambiguity). Section [4] 
discusses identification of principal time-scales and derivation of the limiting models. Section 
[5] reviews general averaging methods, and Section [6] gives additional examples. We believe 
that these methods provide tools for the systematic reduction of highly complex models. 



Further evidence for that claim is provided in Kang (2009) in which the methods are applied 



to obtain a three time-scale reduction of a model of the heat shock response in E. coli given 



by Srivastava, Peterson, and Bentley (2001). 



1.1 Terminology 

This paper relies on work in both the stochastic processes and the chemical physics and 
biochemical literature. Since the two communities use different terminology, we offer a brief 
translation table. 



Chemistry 



Probability 



propensity 

master equation 

Langevin approximation 

Van Kampen approximation 

quasi steady state/partial equilibrium analysis 



intensity 
forward equation 
diffusion approximation 
central limit theorem 
averaging 



The terminology in the last line is less settled on both sides, and the methods we will 
discuss in Section [5] may not yield "averages" at all, although when they don't they still 
correspond well to the quasi-steady state assumption in the chemical literature. 
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2 Equations for the system state 

The standard notation for a chemical reaction 

A + B -± C 

is interpreted as "a molecule of A combines with a molecule of B to give a molecule of C." 

A + B ^ C 

means that the reaction can go in either direction, that is, in addition to the previous reaction, 
a molecule of C can dissociate into a molecule of A and a molecule of B. We consider a 
network of reactions involving s chemical species, Si, . . . , S SQ , and r chemical reactions 

so so 

UikSi ~^ Yl u ik S i' k = l,...,r , 

1=1 1=1 

where the v ik and u' ik are nonnegative integers. If the fcth reaction occurs, then for i = 
1, . . . , s , v ik molecules of Si are consumed and v' ik molecules are produced. We write re- 
versible reactions as two separate reactions. 

Let X(t) G N s ° be the vector whose components give the numbers of molecules of each 
species in the system at time t. Let v k be the vector with components vi k and v' k the vector 
with components v' ik . If the fcth reaction occurs at time t, then the state satisfies 

X(t) = X(t-) + u' k - u k . 

If Rk(t) is the number of times that the kth reaction occurs by time t, then 

X(t) = X(0) + Y,Rk(tM-"k)=X(0) + (v' -u)R{t), 
k 

where v' is the Sq x ro-matrix with columns given by the is' k , v is the matrix with columns 
given by the I/*,, and R(t) G W° is the vector with components Rk{t). 
Modeling X as a continuous time Markov chain, we can write 

R k (t) = Y k ([ X k (X(s))ds), (2.1) 
Jo 
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where the Yk are independent unit Poisson processes and Afc(x) is the rate at which the kth 
reaction occurs if the chain is in state x, that is, A&(X(t)) gives the intensity (propensity in 
the chemical literature) for the kth reaction. Then X is the solution of 

X(t) = X(0) + Vn(/ X k (X(s))ds)(u' k -u k ). (2.2) 
k Jo 

Define Ck = v k — Vk- The generator of the process has the form 

Mf(x) = J2^(x)(f(x + Ck)-f(x)). 



Assuming that the solution of (2.2) exists for all time, that is, X jumps only finitely often 
in a finite time interval, 

f(X(t)) - f(X(0)) - f*f{X{8))d8 (2.3) 

Jo 

is at least a local martingale for all functions on the state space of the process X. 



If (2.3) is a martingale, then its expectation is zero and 

ft 



J2f(x)p(x,t) = Y,f(x)p(x,o) + [ 

Jo 



Bf{x)p{x,s)ds, (2.4) 



where p(x,t) = P{X(t) = x}. Taking f(x) = l{ y y(x), (2.4) gives the Kolmogorov forward 
equations (or master equation in the chemical literature) 

p(y, t) = XI Xk ( y ~ ^ p ( y ~ *) _ x k(y)p(y, t). (2.5) 



The stochastic equation (2.2), the martingales (2.3), and the forward equation (2.5) 
provide three different ways of specifying the same model. This paper focuses primarily on 
the stochastic equation which seems to be the simplest approach to identifying and analyzing 
the rescaled families of models that we will introduce. 

In what follows, we will focus on reactions that are at most binary (that is, consume at 
most two molecules), so \k(x) must have one of the following forms: 

Afc Reaction v k 

K' k ->• stuff 

K k Xi Si — > stuff ej 

K k V~ 1 x i (x i — 1) 2Si — > stuff 2ej 

K , k V~ 1 XiXj Si + Sj — > stuff ti + ej 

Here V denotes some measure of the volume of the system, and the form of the rates reflects 
the fact that the rate of a binary reaction in a well-stirred system should vary inversely with 
the volume of the system. Note that if Qk < 0, then Xk(x) must have x\ as a factor. Higher 
order reactions can be included at the cost of more complicated expressions for the A^. 
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Our intent is to embed the model of primary interest X into a family of models X N 
indexed by a large parameter N . The model X corresponds to a particular value of the 
parameter N = N , that is X = X N °. 

For each species i, let ctj > and define the normalized abundance (or simply the 
abundance) for the iVth model by 

Z»{t) = N-^X t N (t). 

Note that the abundance may be the species number (a, = 0), the species concentration, or 
something else. The exponent ai should be selected so that Zf = 0(1). To be precise, we 
want {Zf(t)} to be stochastically bounded, that is, for each e > 0, there exists K et < oo 
such that 

MP{fm V Z»(s)<K e , t }>l-e. 

N s<t 

In other words, we want ctj to be "large enough." On the other hand, we do not want cej to 
be so large that Zf converges to zero as N — > oo. For example, the existence of 8 e such that 

inf P{inf Zf (s) >5 et }>l-e 

N l s<t 1 V ' - ' J _ 

would suffice; however, there are natural situations in which ctj = and Zf is occasionally 
or even frequently zero, so this requirement would in general be too restrictive. For the 
moment, we just keep in mind that ctj cannot be "too big." 

The rate constants may also vary over several orders of magnitude, so we define K k by 
setting n' k = n k NQ k for unary reactions and k^V' 1 = n k NQ k for binary reactions. The (3 k 
should be selected so that the n k are of order one, although we again avoid being too precise 
regarding the meaning of "order one." For a unary reaction, the intensity for the model of 
primary interest becomes 

and for binary reactions, 

^k^ XiXj — Nq K, k ZiZj = N^ k k HfcZiZj 

and 

^V^Xiixi - 1) = NS k+2ai K kZi ( Zi - N^) = N^ +Uk - a K k z t (z t - N^). (2.6) 
The Nth model in the scaled family is given by the system 

Zf(t) = Zf (0) + J2 N ~ a 'M f N^ + ^X k {Z N {s))ds){v' ik - v ik ). 
k Jo 

For binary reactions of the form 2Si — > stuff with a, > 0, X k (z) = K k Zi(z,i — N~ ai ) depends 
on N, but to simplify notation we still write \ k rather than A^. 

Let A N = diag(A r ~ Ql , . . . , N~ aa o), p k = f3 k + v k ■ a, and ( k = v' k — v k . The generator for 
Z^is 

M N f(z) = J2N pk \ k (z)(f(z + A N ( k ) - f{z)). 
k 
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Even after the (3 k and a» are selected, we still have the choice of time-scale on which to 
study the model, that is, we can consider 

Z?>\t) = Z»{tW) = Z» (0) + J2 N ~ ai M f N^ + ^ a \ k {Z N >\s))ds){v' lk - v ik ) (2.7) 

u JO 



for any 7 G 1. Different choices of 7 may give interesting approximations for different 
subsets of species. To identify that approximation, note that if limAr^oo Z^ a = Z] and N Q 
is "large" , then we should have 

X^^X^^^N^ZUtN^). 



In what we will call the classical scaling (see, for example, Kurtz (1972, 1977/78)) N 
has the interpretation of volume times Avogadro's number and aj = 1, for all i, so Zj~° is the 
concentration of Si. Taking (3 k = for a unary reaction and (3 k = —1 for a binary reaction, 
the intensities are all of the form NX k (z), and hence taking 7 = 0, Z N = Z N '° converges to 
the solution of 



u JO 



; ik -"ik), (2.8) 

where z Vh = Yli z i ik - Note that ( |2.8[ ) is just the usual law of mass action model for the 
network. 



3 Determining the scaling exponents 

For systems with a diversity of scales because of wide variations in species numbers or rate 
constants or both, the challenge is to select the oti and the (3 k in ways that capture this 
variation and produce interesting approximate models. Once the exponents and Nq are 
selected, 

a?(0) = L(0* *(0>J, 

and the family of models to be studied is determined. 
Suppose 

/ \ / \ ^ / 

K l d. K 2 — ' ' ' — K r - 

Then it is reasonable to select the f3i so that 0i > • • • > ro , although it may be natural to 
impose this order separately for unary and binary reactions. (See the "classical" scaling.) 

Typically, we want to select the so that Z^(t) = N~ ai Xf {t) = 0(1), or more precisely, 
assuming linijv->.oo 2^(0) = Zi(0) > 0, for all i, we want to avoid a, (3, and 7 for which 
liniTv-s.oo Zf (iiV 7 ) = 0, for allt > or liniA^oo Zf (tiV 7 ) = 00, for all t > 0. This goal places 
constraints on a, (3, and possibly 7. 

3.1 Species balance 

Consider the reaction system 



S 1 + S 2 

s 3 + s 5 



5*3 + S4 
Dr. 
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Then the equation for is 

Z^{t) = Z£(0) + Ar- a 3y 1 (Ar7+ft+«i+«2 f Kl z^( s )Z^(s)ds) 

Jo 

-N- as Y 2 (N^ 2+a3+a5 [ K 2 Z^ 1 (s)Z^ n (s)ds) . 

Jo 

Assuming that Zf ' 7 = 0(1) for i ^ 3 and Z*(0) = 0(1), Zf' 7 = 0(1) if 

(Pi + ai + a 2 + 7) V (/3 2 + a 3 + a 5 + 7) < a 3 
(the power of N outside the Poisson processes dominates the power inside) or if 

Pi + ai + a 2 = (5 2 + a 3 + a 5 . (3.1) 



Assuming (3.1), if Z^ > KlZl '^fe'^ , the rate of consumption of S3 exceeds the rate 

of production, and if the inequality is reversed, the rate of production exceeds the rate of 
consumption ensuring that Z% neither explodes nor is driven to zero. 

In general, let Tf = {k : v' ik > u ik }, that is, Tf gives the set of reactions that result in 
an increase in the ith species, and let T~ = {k : v' ik < z/^}. Then for each i, we want either 

max(/3 fc + Uk ■ a) = max(/3 fc + u k ■ a). (3.2) 
fcerr fc e r+ 



or 



max (p k + u k ■ a) + 7 < a { . (3.3) 
fcer+urr 



We will refer to (3.2) as the balance equation for species i and to (3.3) as a time-scale 
constraint since it is equivalent to 

7 < a, - max (P k + v k -a). 
fcer+urr 

The requirement that either a species be balanced or the time-scale constraint be satisfied 
will be called the species balance condition. 

Equation (3.2) is the requirement that the maximum rate at which a species is produced 
is of the same order of magnitude as the rate at which it is consumed. Since consumption 
rates are proportional to the normalized species state Zi, Z{ should remain 0(1) provided 
the same is true for the other Z 3 - even if the normalized reaction numbers blow up. If (3.2) 
fails to hold, then (3.3) ensures that Zi(t) = 0(1), again provided the other Zj remain 0(1). 

Note that if ( ik ^ 0, then 

7 = lik = on - (P k + v k ■ a) (3.4) 
is in some sense the natural time-scale for the normalized reaction number 

N-^R^(t) = N^Y k (N^ k+Uk - k [ Xk(Z N ^(s))ds). 

Jo 



Then, regardless of whether (3.2) or (3.3) holds, 



Ti 



min j ik = oti- max ((3 k + v k ■ a) 



(3.5) 



is the natural time-scale for species Si. With reference to (2.7), if 7 < 7*, we expect Z J Ar ' 7 (t) 
to converge to hin^^oo Z^(0). If 7 = 7$ and ctj > 0, then we expect 

lim Zf'^(t) = Km (Zf (0) + £ f \ k (Z N ™(s))ds(u' ik - v ik % 



where 



T ifi = {I : & + vi ■ a = max (/3 fe + z/ fc ■ a)} 

fcer+urr 



and each integral on the right side is nonconstant but well behaved. If aj = 0, we expect 
lim Z***(t) = lim {Z? (0) + £ n( t \ k (Z N ^(s))ds)(v' ik - u lk )). 



It is important to notice that we associate "time-scales" with species (and as we will see 
below, with collections of species) and that one reaction may determine different time-scales 
associated with different species. 



3.2 Collective species balance 

The species balance condition, however, does not by itself ensure that the normalized species 
numbers are asymptotically all 0(1). There may also be subsets of species such that the 
collective rate of production is of a different order of magnitude than the collective rate of 
consumption. Consider the following simple network: 

If < 04 < (3i < /3 2 = /?3 and a.\ = a 2 = 0, then 

Zf (*) = Zf(0) + Y^NN) + Y 3 (k 3 N^ f Z 2 N (s)ds) - Y 2 (k 2 N^ [ Z^(s)ds) 

Jo Jo 

Z»(t) = Z?(0) + Y 2 (k 2 N^ f Z?(s)ds)-Y 3 (K 3 N^ [ Z^(s)ds) (3.6) 

Jo Jo 

-Y 4 {k 4 N^ 4 [ Z^(s)ds) . 
Jo 

Since {3 2 = (3 3 V {3\ and {3 2 = (3 3 V /3 4 , the species balance condition is satisfied for all species, 
but noting that 

Zf (t) + Z?(t) = Z? (0) + Z?{0) + Y^kxNN) - Y a (k a N^ f Z^(s)ds), 

Jo 
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the species numbers still go to infinity as N — > oo. This example suggests the need to 
consider linear combinations of species. These linear combinations may, in fact, play the role 
of "virtual" species or auxiliary variables needed in the specification of the reduced models 



(cf. |Cao, Gillespie, and Petzold] 02005) and |E, Liu, and Vanden-Eijnden| p005| |2007[ )). 
To simplify notation, define 

Pk = (3k + Vk-a, 

so the scaled model satisfies 

Z N "(t) = Z N ^{0) + A N V Y k {N^ +u "- a+ ^ f \ k {Z N ~<{s))ds)( k 

= Z N ^(0) + A N J2 Yk(N pk+1 [ \ k (Z N "(s))ds)( k , 
k Jo 

where A^v is the diagonal matrix with entries N~ ai . 

Definition 3.1 For 9 G [0, oo) s °, define Tg={k:9-( k >0} and Tg = {k : 9 ■ ( k < 0}. 
Then, noting that 



so 



so 



i=l 



i=l 



9 T A~ N 1 Z N ^{t) = e T A N 1 Z N '\O) + J2(0-Ck)Y k (N^' f \ k {Z N ^{s))ds) 

k Jo 

= e T A- N 1 z N '-r(o) + • Ck)R^(t) - £ |(* • Ck)\R^(t). 



To avoid some kind of degeneracy in the limit, either the positive and negative sums must 
cancel, or they must grow no faster than N ai for some % with 0i > 0. Consequently, we extend 
the species balance condition to linear combinations of species. For each 9 G [0,oo) s °, the 
following condition must hold. 



Condition 3.2 



max(/3 fc + v k ■ a) = max(/3 fc + v k ■ a) 



fcer fl 



ker^ 



or 



7 < = max at - max (f3 k + v k ■ a). 
v.Qi>$ fcer+ur" 



(3.7) 
(3.8) 



Of course, if 9i > for only a single species, then this requirement is just the species 



balance condition, so Condition 3.2 includes that condition. Again, we will refer to (3.7) 



as the balance equation for the linear combination 9 ■ X = 9iX{. In the special case of 
9 = ej, the vector with zth component 1 and other components 0, we say that Xi is balanced 
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or that the species Si is balanced. If (3.7) fails for 9, we say that 9 ■ X is unbalanced. The 

(3.9) 



inequalities given by (3.8) are again called time-scale constraints as they imply 

7< 



mm 70. 
e x unbalanced 



For example, consider the network 



«3 



and assume that = KkN^ k , where {3± = fl 2 > p3- For S 2 to be balanced, we must have 
P2 + oil = + 0L2 and for Si to be balanced, we must have 



/?i V % + a 2 ) = (3 2 + a 



i- 



Let ai = and a 2 = (3 2 — (3 3 so S x and S 2 are balanced. For 9 = (1, 1), = {1}, and 
T e = 0. Consequently, (3.7) fails, so we require 



7 < «i V a 2 - f3x = -/3 3 - 



(3.10) 



There are two time-scales of interest in this model, 7 = — f3\, the natural time-scale of Si 
and 7 = — (3 3 , the natural time-scale of S 2 . The system of equations is 

-f 



Zf' 7 (t) 



zf^ + yx^iv^ 1 *)-^^^ 2 / zf' 7 ( s )rf s ) 

Jo 

+Y 3 (k 3 N^ 3 



Z% {$) + N- a *Y 2 {K 2 N^ [ Z^\s)ds) 

Jo 

-N~ a2 Y 3 (K 3 N^ +a2 I Z 2 a (s). 

Jo 

For 7 = — /?!, since /?i = /3 2 = /3 3 + a 2 , the limit of Z^' 7 satisfies 

Z x {t) = Z 1 (0) + Y 1 ( Kl t)-Y 2 (K 2 [ Z 1 (s)ds) + Y 3 (K 3 [ Z 2 {s)) 

Jo Jo 



= Z 1 {Q) + Y l {K 1 t)-Y 2 {K 2 Z l (s)ds) + Y 3 (K 3 Z 2 (0)t) 

Jo 

Z 2 (t) = Z 2 (0). 

For 7 = —f3 3 , if we divide the equation for Z^ a by N a2 = iV^ 1-/33 , we see that 

= lim N~ a2 Z?'\t) (3.1i; 

N— >oo 



lim N~ a2 Z^((}) + N-^Y^N^H) - N~ a2 Y 2 (K 2 N^ +l32 [ Z^ 1 (s)ds) 

+N- a2 Y 3 (K 3 N^ 3+a2 [ Zf' 7 (s) 

Jo 

lim I Kit + k 3 Z 2 ' 1 {s)ds — k 2 j Z^ )1 {s)ds ) 
V Jo Jo J 
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and Z 2 a converges to 

Z 2 {t) = Z 2 {0) + Kit. 

With reference to ( |3.10[ ), if 7 > -/3 3 , then Z 2 ^(t) —7- 00, for each t > 0, demonstrating the 
significance of the time-scale constraints. 

For 7 = — fa, Z^ 1 fluctuates rapidly and does not converge in a functional sense. Its 
behavior is captured, at least to some extent, by its occupation measure 



V^(Cx[0,t}) 



l c (Z N "(s))ds. 



Applying the generator to functions of Z\ and using the fact that {3\ — (3 3 = (3 2 — (3 3 = a 2 , 

B N >y{z u z 2 ) = N a *C z J{ Zl ), where 



C z J{z x ) = («! + K 3 z 2 )(f( Zl + 1) - f{z x )) + K 2Zl {f{z x - 1) - f[z x )). 



Then 



f{Z?"{t)) - /(^ v ' 7 (0)) - N°* / C^, 7(s) M)VV v ' 7 (^i x ds) 

Jm[o,t] 2 



is a martingale, and dividing by N a2 and passing to the limit, it is not difficult to see that 
V-^'" 1 converges to a measure satisfying 



/ 



Nx[0,i] 



C^ w /(^i)Vi(d^i x ds) = 0. 



(See Section[5j) Writing Vx{dz\ xds) = v s (dzi)ds, it follows that v s is the Poisson distribution 
with mean Kl+K3Z2 W _ We will refer to v s as the conditional-equilibrium or local- averaging 



distribution. 



3.3 Auxiliary variables 



While (3.5) gives the natural time-scale for individual species, it is clear from examples 



considered by E, Liu, and Vanden-Eijnden (2005), that the species time-scales may not be 



the only time-scales of interest. For example, they consider the network 



t\> -1 t\> o Aj rr 

^ Q ^ Q °^ Q 

1 T O4 
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with k[, k' 2 , k' 5 , k' 6 >> k' 3 , k' 4 . The scaled model is given by 

Z?(t) = Z^(0) + N~ ai Y 2 (K 2 N^ +a2 [ Z 2 N (s)ds) - iV- ai r 1 (/t 1 iV /3l+ai / Z?(s)ds) 

Jo Jo 

Z?(t) = Z?(0) + N- aa Y 1 ( Kl N f>1+ai [ Z*{s)ds)-N~ a2 Y 2 {K 2 N P2+a2 [ Z* \s)ds) 

Jo Jo 

+N~ a2 Y 4 (n 4 N l34+a3 [ Z 3 N (s)ds) - N- a2 Y 3 (K 3 N k+a2 [ Z* \s)ds) 
Jo Jo 

Z*(t) = Z^(0) + N- a3 Y 6 (K 6 N^ e+a4 [ Z^(s)ds) - N~ a3 Y 5 (K 5 N^ +a3 [ Z^(s)ds) 

Jo Jo 

+N~ a3 Y 3 (K 3 N P3+a2 [ Z^(s)ds) - N~ a3 Y 4 (K 4 N^ +a3 [ Z 3 N (s)ds) 
Jo Jo 

Zf(f) = Z?{0) + N- a4 Y 5 {K 5 N^+«3 f Z* \s)ds) - N- a4 Y 6 {K 6 N^ +a4 [ Z^(s)ds). 

Jo Jo 

Assume that Pi — P 2 > P5 — Pq > p 3 > P 4 - Then if we look for a scaling under which all 
6 • X are balanced, oti = a 2 , a 3 = a 4 , and a 2 + (3 3 = a 3 + (3 4 , so a 3 = a 2 + (3 3 — p 4 . For 
definiteness, take 0:1 = 0:2 = 0. 

The natural time-scale for Si and S 2 is —pi, and the natural time-scale for S 3 and S 4 is 
—p5, but on either of these time-scales Z\ + Z 2 and Z 3 + Z 4 are constant. In particular, 

U^(t) = Z^\t) + Z^'\t) = Z? (0) + Z? (0) + Y 4 (n 4 N^ +a3 f Z 3 N '\s)ds) 

Jo 

-Y 3 (k 3 N^ 3 f 'z^(s)ds) 
Jo 

U 2 '\t) = Z^(t) + Zf' 7 (t) = Z?(0) + Z? (0) - N~ a3 Y 4 ( K4 N^ +a3 [ Z 3 Nr/ (s)ds) 

Jo 

+N~ a3 Y 3 (K 3 N^ 3 [ Z^ n {s)ds). 
Jo 

For 71 = 72 = -pi = -p 2 , {Z^ ai ,Z 2 ' lr ) converges to 

Z~{\t) = Zi(0)+Y 2 (k 2 f Z^ 1 (s)ds)-Y 1 (K l f Zj 1 (s)ds) 

Jo Jo 

Z?(t) = Z 2 (0) +Yi(ki fzl\s)ds)- Y 2 (k 2 f Z?{s)ds), 

Jo Jo 

and for 73 = 74 = ~p 5 = -/3 6 , 

Zf{t) = Z 3 (0) + k 6 f Zj 3 (s)ds-K 5 f Zf{s)ds 

Jo Jo 

Zf{t) = Z 4 (0) + k 5 [ Z^ 3 (s)ds-K 6 [ Zf(s)ds. 

Jo Jo 
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Taking 712 = — 03 = —(03 + ^4) and dividing the equation for Z^' 712 by /?3 , we see 
that 



k 5 [ Z^ lu (s)ds- k 6 [ Z? r(12 (s)ds ->• (3.12) 
Jo Jo 



and hence ^ ^ 

Z^ 12 (s)ds — f U^ 12 (s)ds ->0. (3.13) 

^5 + ^6 Jo 







Similarly, dividing the equation for Z 1 ' 112 by N 132 133 , 



* Zf 7l2 (s)ds ^— f C/f ' 7l2 (s)ds 0. 

Since U^ r/12 converges to (0) uniformly on bounded time intervals, [7^' 712 converges to 
the solution of 

Ur(t) = tMO) + Y 4 (J^U 2 (0)t) - Y 3 (^^- f U{ 2 {s)ds). 

K 5 + K 6 Ki + K 2 Jo 



Finally, taking 734 = —pi, as in (3.13), 

f Z^ 4 (s)ds ^— f U^ 34 (s)ds 0, 

JO K 5 + K 6 JO 

and dividing the equation for JJi'" 123 by jV^ 3 ™^ 4 , 

f z^ s4 { S )ds - — r ^' 734 ( S )d S -> 0. 

Jo k 3 Jo 

Consequently, even on this faster time-scale, U2'™ converges to £72 (0) uniformly on bounded 
time intervals. 

3.4 Checking the balance conditions 



Condition 3.2 only depends on the support of 9, supp(6>) = {% : 9{ 7^ 0}, and on the signs 
of 9 ■ (k, so the condition needs to be checked for only finitely many 9. For k G {1, . . . ,r }, 
define 

A+ = {9 G [0, oo) so : 9-C k > 0}, A, = {9 G [0, oo) s ° : 9-( k < 0}, A° = {9 G [0, 00)*' : 0< fc = 
and for disjoint r_, T + , T satisfying r_ U T + U T = {1, • • • , r }, define 

A r _,r + ,r = (n te r_A fc ) n (n fcer+ A+) n (n fce r A°). 
The following lemma is immediate. 
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Lemma 3.3 Fix 7. Condition 3.2 holds for all 9 £ [0, oo) s ° provided 



max(/3 fc + u k ■ a) = max(/3 fc + v k ■ a) 



(3.14) 



or 



7 < min max oti — max (B k + u k ■ a) 

se^r_T + ,r i.di>o fcer + ur_ 



(3.15) 



for all partitions {r_,r + ,ro} for which Ar_,r +I r 7^ 



Checking the conditions of Lemma 3.3 could still be a formidable task. The next lemmas 



significantly reduce the effort required. Observe that for 9 1 ,9 2 £ [0,oo) s ° and C\,C2 > 0, 
k £ r+ i6)1+C2(9 2 implies k £ rji U T+ 2 and similarly for r~ el+C202 , so 



max pfc < 



< max pfc V max p k 



(3.16) 



and 



max pfc < max p k V max p^. 



fcer cl ei +C2 e 



fcer 



(3.17) 



Let G be a directed graph in which the nodes are identified with the species and a 
directed edge is drawn from Si to Sj if there is a reaction that consumes Si and produces 
Sj. A subgraph Go C G is strongly connected if and only if for each pair Si, Sj £ Go, there 
is a directed path in G beginning at Si and ending at Sj. Single nodes are understood to 
form strongly connected subgraphs. Recall that G has a unique decomposition G 
into maximal strongly connected subgraphs. 



The following lemma may significantly reduce the work needed to verify Condition 3.2 



Lemma 3.4 Let 9 £ [0, oo) s °, and fix 7. Write 



(3.18) 



where supp(^) C Gjfor some maximal strongly connected subgraph Gj and Gj 7^ G,, /or 
z ^ j. If Condition 3.2 holds for each Qi , then it holds for 6. More specifically, if the balance 



equation (3.1) holds for each 6 3 , then the balance equation holds for 6, and if (3.8) holds for 



each 3 , then (3.8) holds for 6. 



Consequently, if Condition 3.2 holds for each 6 £ [0, oo) so with support in some strongly 
connected subgraph, then Condition 3.2 holds for all 9 £ [0,oo) so ; if (3.1) holds for each 
9 £ [0, oo) s ° with support in some strongly connected subgraph, then (3.1) holds for all 
9 £ [0, oo) s °; and if (3.8) holds for each 9 £ [0, oo) s ° with support in some strongly connected 
subgraph, then (3.8) holds for all 9 £ [0, oo) s °. 



Proof. Assume that Condition 3^ holds for each 9- \ j 
7^ 0. Select l\ £ satisfying 



Ph 



max p k . 

fcer+ 



1, . . . , m. First, assume that 

(3.19) 
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Since Tq C UjTj~, there exists ji such that l\ G ri i; and using (3.19), we have 

maxp/c = p/ x < max pk- (3.20) 



We have three possible cases. First, if max, pr + pk 7^ max fcpr - pk, then by (3.8), there 

eh eh 



exists i\ G supp(6 n ) such that 



7+ max Pk<oc h , (3.21) 



fcer + ur _ 
en en 



and by (|3_20j), 

7 + maxpfc < ojj, < max «j. (3.22) 
feer+ iesupp(e) 



Second, if max fcgr + pt = max fcgr - pk < max fcgr - pk, then by (3.20), we obtain 

eh eh « 



maxpfc < max pk = max pk < maxpfc. (3.23) 

Finally, if 



fcG r+ ker+ n ker^ k& g 



max pk = max pk > max pk (3.24) 



feer+. fcer" fcer" 

en en 



we select l 2 in r g3l with p/ 2 = max^p- p fc . The fact that pi 2 > max fcgr - p k ensures the 



en 



existence of j% such that l 2 G r+ 2 . Then we have 



max pk = max pk = pi 2 < max pk- (3.25) 



fcer + fcer - . ker^ 

en en en 

We recursively select l n and j n with /„ G r^ n such that 

max p k = max p fc = p in < max p fc 
fcer+. feer - . fcer+ 

until we find Z n for which this is no longer possible. Since the Gj are maximal strongly 
connected subgraphs, there is no possibility that the same Qi is selected more than once. 
Thus, the process will terminate for some n and when it does max, pr + pk 7^ max, er - pk 

ein gin 

and 



7 + maxpfc < 7 + max pk < max a i — a i- (3.26) 

fcer^ fc€ 

Consequently, we always have either 



fcer+ fcer+ n jesupp(0J") iesupp(e) 



7 + maxpfc< max ai (3.27) 
fcer + iesupp(e) 

or 

maxpfc < maxpfc. (3.28) 
feer+ fcerv 
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If T e 7^ 0, interchanging — and +, we see that either 



or 



7 + max pk < max on 
feG r- iesupp(e) 



max pk < max pk . 
fcer- fcer+ 



(3.29) 
(3.30) 



Assume that both Tg and T e are nonempty. If both (3.28) and (3.30) hold, then (3.7) 
is satisfied. If (3.27) and (3.29) hold, then taking the maximum of the left and right sides, 



(3.8) holds. If (3.27) and (3.30) hold, then (3.8) holds and similarly for (3.28) and (3.29). 



If (3.7) holds for all 9 3 , then the first and third cases above cannot hold so (3.23) must 
hold giving (3.28) and by the same argument (3.30). Consequently, (3.7) must hold for 9. If 
(3.8) holds for all 9 J , then the first case above holds giving (3.27) and by the same argument 
(3~29l, so ([3781 must hold for 9. 



If T~q = and Tg ^ 0, then (3.29) must hold and (3.8) holds for 9 and similarly with the 
+ and — interchanged. 

If both and Tq are empty, then (3.7) holds (— oo = — oo). In particular, 9 ■ (k = for 
alia- ' □ 



The remaining lemmas in this section may be useful in verifying Condition |3.2| for the 
cases that remain, that is, for 9 with support in some strongly connected subgraph. 



Lemma 3.5 Fix 7 G R, and suppose (3.8) holds for 9 1 ,. . . ,9 m E [0, oo) s °. Then for Cj > 0, 
j = 1, . . . , m, (Q holds for 9 = Y^=i Cj0 j . 

Proof. Since 9 ■ (k > implies Cj9 J ■ (k > for some j and 9 ■ (k < implies Cj9 J ■ (k < for 
some j, 

max pk < max max pf. 
feer+ur" 1 <J<"»jfe6r+.ur; 



and there exists j such that 



7 < max oti — max p& < max ctj — max p^. 
i-.e{>0 ker tj ur 9~ 3 iA>0 feer e ur V 



□ 



Lemma 3.6 For9\9 2 e [0, oo) s ° ; suppose that 



max pk = max pk > max pk- 
fcer- fcer+ fc er+ur- 2 



(3.31) 



Then for c u c 2 > 0, (3. 1) holds for Cl 9 l + c 2 9 2 



Proof. If I G and p\ = max fcgr + pk, then by (3.31), / ^ T d2 . Consequently, / G r^ 01 
and by (3.16), we must have 



fcer 



max pk = max pk- 
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By the same argument, 



max 



Pk 



and it follows that (3.7) holds for c\9 l + c 2 9 2 . 



max pk, 



□ 



Lemma 3.7 Fixj, and suppose that (3.1) holds for 9 1 and (3.8) for 9 2 . Then for C\,c 2 > 0, 
Condition\3.2\ holds for Cy9 x + c 2 9 2 . 



Proof. If 



max pk = max pk > max pk V max pk, 



ker 



fcer: 



then Lemma 3.6 implies C\9 X + c 2 # 2 satisfies (3.7), so assume that 



max pk = max pk < max p^ V max p^. 
fcer «"i fcer ^"i fcer « _ 2 fc6r+ 



(3.32) 



(3.33) 



Then 



and 



so 



ker 



max pk < max p& V max pk < max pk V max p^ 



lS 1 +c 2 9 2 



fcer 



fcer 



fcer 4 



max pfc < max pk V max p^. < max pk V max pk, 



fcer 



fcer 



ker 



max pfc V max pk < max pk V max p^, 

c 1 9 1 +c 2 e 2 



ker: 



feer -1 



and since supp(ci6 |1 + c 2 # 2 ) D supp(# 2 ), (3.8) for 9 2 implies (3.8) for c\9 x + c 2 9 



□ 



If Condition 3J2 holds for 9 1 and 8 2 and c\,c 2 > 0, then the previous lemmas imply 

max pk = max pk = max pk = max p^. (3.34) 



Condition 3.2 holds for ci^ 1 + c 2 9 2 except in one possible situation, that is, 



fcer fl "i 



ker 



fcer 



Since the species balance condition does not imply Condition 3.2 for 9 = (1, 1) for the 
system (Zf ,Z^) given by (3.6), some additional condition must be required to be able to 



conclude Condition 3.2 holds for C\9 X + c 2 9 2 when (3.34) holds. The following lemmas give 
such conditions. 



Lemma 3.8 Fix 7 G E, and suppose that Condition 3.2 holds for 9 1 ,9 2 G [0,oo) s °. If 
Tji fl = or r fll n r+ 2 = and ci, c 2 > ; £/ien Condition 3.2 holds for c\9 l + c 2 9 2 . 
If Condition 3.1 holds for 9\ and 9 2 , Tti fl T e2 - 



or r~! n r 



Condition 3.1 holds for C\9 X + c 2 9 



), and C\,c 2 > 0, then 



Remark 3.9 If no reaction that consumes a species in the support of 9 1 produces a species 
in the support of 9 2 , then F7i D rj~ 2 = 0. That condition is, of course, equivalent to the 
requirement that a reaction that produces a species in the support of 9 2 does not consume a 
species in the support of 9 l . 
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Proof. As noted, the previous lemmas cover all possible situations except in the case that 
(|3.34|) holds. Suppose T el n T+ = 0. If 6 1 ■ ( k < 0, then 6 2 ■ ( k < and {c x B x + c 2 6 2 ) ■ ( k < 0, 



and if (c^ 1 + c 2 # 2 ) • Cfc < 0, then either 6 1 ■ ( k < or 9 2 ■ ( k < 0, so 

max p k < max p k < max p k V max p^. 



fcer «"i 



C 1 B + c 2 6 



fcer; 



(3.35) 



Similarly, noting that # 2 ■ ( k > implies 9 l ■ Q k > 0, 



max p k < max p/c < max p k V max p k . 



(3.36) 



e 2 



But (3.34) implies equality holds throughout (3.35) and (3.36) and (3.7) holds for c\Q x + C2O 2 . 
□ 



Lemma 3.10 Suppose (3.1) holds for 9 1 and 6 2 and for 9 1 - e ^6 2 for all k G (r+ n T 92 ) U 
(r-j n T+ ). (Note that -f^ > 0.) Then (Q holds for c x d x + c 2 9 2 for all c u c 2 > 0. 



Proof. By Lemma 3.6, we can restrict our attention to the case (3.34), and it is enough to 
consider 6 1 + c6 2 for c > 0. Note that for c sufficiently small, Tq 1+cQ2 D and T Q1+ce2 
and 



max p k = max p k = max p k = max p k . 

ei+ce2 fcer e i keT gi +c0 2 fcer s i 



fcer « + l 



Let 



Co = inf{c : max p k 7^ max p k or max p k 7^ max p^}, 
fcer + , , fcer+ fcer~ „ fcer~ 



and note that for < c < cq, (3.7) holds. If Co < 00, then cq 



by the assumptions of the lemma 



max p k 



max < max p^ = max p k . 



fcer" 
ei+c 



a-' 



> for some k, and 



(3.37) 



But (3.37) can hold only if there exists / + G such that pi+ = max fcgr + p k and Cq = — ^ 

-. Then, for c > Cq, / + G r fll 



and / G r 0i such that p«- = max fcgr - p k and Co — —qt^ 



and / G , and the lemma follows. 



e 1 +c6 2 
□ 



4 Derivation of limiting models 

As can be seen from the examples, derivation of the limiting models can frequently be carried 
out by straightforward analysis of the stochastic equations. The results of this section may 
be harder to apply than direct analysis, but they should give added confidence that the 
limits hold in great generality for complex models. 
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We assume throughout this section that limjv^oo Z^' 7 (0) exists and is positive for all i. 

If 

7 = ri = min7i = min(a i - max (f3 k + v k ■ a)), (4.1) 
* * fcer+ur- 

then \ivn N _ ¥00 Z Na exists, at least on some interval [0, Too) with t^ > 0, and is easy to 
calculate since on any time interval over which sup t<T \Z N,1 {t)\ < oo, each term 

N~ ai Y k ( [* W + **\ k {Z N «(s))ds) 



either converges to zero (if oti > 7 + p k ), is dependent on N only through Z Na (if ajj 
7 + Pk = 0), or is asymptotic to 



\ k (Z N >\s))ds 



(if oti = 7 + p k > 0), since 



lim sup \N- ai Y k (N ai u) - u\ = 0, m > 0. 

The caveat regarding the interval [0,Too) reflects the fact that we have not ruled out "reac- 
tion" networks of the form 2 Si — > 3 Si, Si — >■ which would be modeled by 

X 1 {t) = X 1 {0)+Y 1 { Kl I X x {s){X x {s)-l)ds)-Y 2 {K 2 [ X 1 (s)ds) 

Jo Jo 

and has positive probability of exploding in finite time, if -X"i(0) > 1. 



Theorem 4.1 Let T] = {k : 7 + p k = a*}. For 77 denned 6y gjp ; Z^ 1 Z ri on [0, rj, 

where if oti > 0, 

Z?(t) = Z t {0)+ J2 I h(Z n (s))ds(u' ik -u ik ), 



if ol { = 0, 



and 



Z?{t) = Z t {0) + Y ^ f h(Z^(s))ds)(p' ik 



Too = lim t c = inf{t : sup |Z ri (s)| > c}. 



s<t 



Remark 4.2 By Z N,ri =>■ Z ri on [0, Too), iue mean that there exist r N,n and r" such that 
(Z N ' ri (- Ar N ' n ),r N ' n ) => (Z ri (- Ar"),r") and lim„^oc T n = . 

Proof. Let tat jC = inf{t : sup s< jZ Ar,ri (s)| > c}. The relative compactness of {Z N ' Vl (■ A tat jC )} 
follows from the uniform boundedness of X k (Z N ' ri (■ A Tjv, c ))- Then (Z Ar,ri (- A tat >c ), tat jC ) =>■ 
(Z Tl (- A t c ), t c ) at least for all but countably many c. □ 
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Note that 70 > minj^oTi, so 77 = mine 6 [o )0 o) s o 7e, and Condition 3.2 always holds for 



7 = 77. If Condition 3.2 holds for some 7 > 77, then the balance equality (3.7) must hold 
for all 9 G [0, 00) s ° with 70 = 77. 
Let «e = maxj : 0. >o dj, and define 



8=1 

As noted earlier, the natural time scale for is 

70 = ct e - max p fc . 
fcer+ur~ 

Let Li be the space spanned by Si = {e^ : T- 1 7^ 0}, and let II! be the projection of IR S ° 
onto Li. Let 

K 2 = {9 G [0, oo) s ° : 9 ■ IliCfc = OV/e G U^ 1 } 

and L2 = spanK2, and let II2 be the projection onto L2. Of course, K 2 contains §2 = {e^ : 
= 0}, but as in the example of Section 3.3, it may be larger. The projections ITi and n 2 



are not necessarily orthogonal, but for any x G M s °, x — n 2 x G L^. 
Lemma 4.3 For each x G W° , x — U2X G Lj. 

Proof. Note that L x = {x G M s ° : • x = 0, Ve^ G S> 2 } and that for e^ G § 2! • n 2 a; = e^- x. 
Consequently, for ej G §2, • (x — n 2 x) = and x — U 2 x G L x . □ 

Lemma 4.4 // G K 2 and 9 ■ Q ^ for some I G T^ 1 , t/ien ojg > ojj. 
Lei 

r2 = min76i = minia^ — max p^}. 

eeK 2 e 6 K 2 L feer+ur- 

T/ien r 2 > ri . 

Proof. For G K 2 , if 9 ■ Q 7^ 0, then 9j > for some j such that e^- ^ Li. Since 



r i < ocj — max pjt < ojj — p/ 
feer+urj 

and T^ 1 = {A; : 07 — pk = 77} = 0, 

n = aci- pi < 07 - p h 



and ae > 07 > at{. 

Let G K2 satisfy 70 = r 2 . Then there exists 



I g r+ u r, c u i:flj>0 rt u rj 



and 0j > such that i G f| U r • and 



r 2 = a fl - # > aj - Pz > 7j > r i- 
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If atj — pi > 7j or jj > ri, then r 2 > r%. If <x,- — pi — 7j = ti, then / e T^ 1 and since 6 1 • Q ^ 0, 
> a, and r 2 > r x . □ 

Unfortunately, while r 2 can naturally be viewed as the second time scale, we cannot 
guarantee a priori that the system will converge to a nondegenerate model on that time 
scale. For example, consider the network 



->• S 1 -> S 2 -> S 3 
5a + S 2 -> Si + S 3 -> 
and assume that the parameters scale so that 

X 1 (t) = X 1 (0) + Y l (K 1 t)-Y 2 (K 2 [ X 1 (s)X 2 (s))ds)-Y 5 (K 5 N- 1 [ X l (s)X 3 (s)ds) 

Jo Jo 

X 2 (t) = X 2 (0) + Y 3 (n 3 t)-Y 2 (K 2 [ X 1 ( 8 )X 2 ( 8 ))d8) 

Jo 

X 3 (t) = XsM+YifaN-ty-YsfaN- 1 I X 1 (s)X 3 (s)ds). 

Jo 



Then Condition 3.7 is satisfied for all 9, r\ = 0, and r 2 = 1. But if K\ > k 3 , Xi(Nt) — > oo 
and X 3 (Nt) ->■ for all t > 0. 

The problem is that even though the balance equations are satisfied for the fast sub- 
network (Xi,X 2 ), the subnetwork is not stable. Consequently, to guarantee convergence 
on the second time scale, we need some additional condition to ensure stability for the fast 
subnetwork so that the influence of the fast components can be averaged in the system on 
the second time scale. 

Of course, with reference to (3.11) and (3.12), it is frequently possible to verify conver- 
gence without any special techniques, but we will outline a more systematic approach. 
We assume the following condition on the scaling. 



Condition 4.5 For each N, AtvIT 2 A 



-i 

N 



n 



2- 



Let L a be the span of {e, : a% = a}. Then Condition 4.5 is equivalent to the requirement 
that II 2 : L Q ->■ h a . 

Define the occupation measure on L x x [0, oo) by 

^ r2 (Cx[0,i])= f 1 C {{I -Yi 2 )Z N ^{s))ds. 
Jo 



Assume that 
in the sense that 



N,r 2 



Vi 



(4.2) 



f{x)V" >ra {dx x ds) 



Li x [0,t] 



f(x)Vi(dx x ds) 



Li x [0,t] 
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for all / G Cfc(Li) and all t > 0. This requirement is essentially an ergodicity assumption on 
the fast subsystem. 

Define = inf{t : \U 2 Z N ' r2 {t)\ > q}. For 9 E K 2 , define T r e 2 = {k:r 2 + p k = a e } and 

h q ,e{y)= sup ^2 \6 ■ ( k \\ k (x + y), 

x£L 2 ,\x\<1 k&v r 2 



if 



and assume that for q > 0, : [0, oo) — > [0, oo) satisfies lim^oo r ^^(r) = oo and 

^eMy^V^idyxds)} (4.3) 

vf] 

is stochastically bounded. In addition, assume 

V \9 ■ Ck\N r * + ^ a ° [ X k (U 2 Z N ^{s) + y)V^ r \dy x ds) -> 0. 

Then for all but countably many q, at least along a subsequence, U 2 Z N,T2 (- A rj^) converges 



in distribution to a process Z r 2 2 {- A r g ), and if p k + r 2 = a$, by Lemma A. 6 



X k (Z N > r2 (s))ds^ / X k (Z r 2 2 (s)+y)V 1 (dyxds). (4.4) 

JLix[0,tAr,] 

Theorem 4.6 Define D a = diag(. . . l{ ai=a y . . .). Under the above assumptions, there exists 
a h 2 -valued process Z r 2 and a random variable > such U 2 Z N,r2 converges in distribution 
to Z r 2 on [0, Too). For 9 e ~K 2 with otg = 0, 

e ■ z?(t) = e ■ z r 2 2 (o) + V (e ■ c*)n( / h(z?(s) + y )v x {d y x ds)) 

fcer r 2 ^LlX[0,t] 

and /or # G K2 with ag > 0, 

9 ■ D a °Z r 2 2 (t) = 9 ■ D a »Z?(Q) D ° a Ck) [ MZ r 2 2 {s) + yW^dy x ds), 

k( _ r r 2 </LlX[0,i] 

forte [O^oo). 

In particular, Z^ ,r2 =>• 9 ■ D ae Z r 2 2 . 

Remark 4.7 The statement of this theorem is somewhat misleading. Given V\, Z 2 2 is 
uniquely determined. However, as we will see in the next section, typically V% depends on 
Z r 2 2 . There we will give conditions under which the sequence of pairs {{V^ r,r2 , Z N,r2 )} is rela- 
tively compact. Then any limit point (Vi, Z 2 2 ) will satisfy the equations given by the present 
theorem, but it will still be necessary to show that the pair is uniquely determined. 



Proof. As for the first time-scale, stopping the process at 

rf = M{t:\U 2 Z N ^{t)\>q) 

ensures that {U 2 Z N,T2 (- Ar^)} is relatively compact, and (4.4) ensures that any limit process 
satisfies the stochastic equations. Uniqueness for the limiting system then follows by the 
smoothness of the X k . □ 
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5 Averaging 



Stochastic averaging methods go back at least to Khas'minskii ( 1966afb ). In this section we 
summarize the approach taken in Kurtz (1992). See that article for additional detail and 
references. 

Recall that Ajv = diag(iV~ ai , . . . , N~ as o) y p k = f3 k + i / k -a, and ( k = v ' k — v k . The generator 
for Z N >° is 

M N f(z) = J2N pk \ k (z)(f(z + A N ( k ) - f(z)). 

k 

Another way of characterizing 7*1 is as the largest 7 (possibly negative) such that limAr^oo N' y Mpff(z) 
exists for each / G C c 2 (M m ) and z G M m Define T 1 ^ = {k : r x + p k = a} and set 
D a = diag(. . . l K=a} . . .). Then 



C Q f(z) 



lim N ri M N f(z) 



h(z)(f(z + A°C») - f{z)) + (E E X ^)D a (k) ■ V/(; 



fcer. 



which is the generator for the limit of the system on the first time scale. The state space for 
the limit process is E = YliLi where Ej = N if aj = and Ej = [0, 00) if Qjj > 0. 

Note that if G T^ 1 , then D a ( k = Tli(k, and by the definition of L 2 , n 2 IIi^fc = 0. 
Consequently, for z G n 2 E and 

E z = {y G Li : y = (I - U 2 )x, Il 2 x = z, x G E}, 

C'f{y) = C„/(z + y) 

defines a generator with state space E z . 
As before, define 

^ r2 (Cx[0,i])= f \ C ((I -U 2 )Z N ^(s))ds, 
Jo 

and observe that 

Mf{t) = f{Z N ' r *{t)) -f(Z N ' r2 {0)) - [ N r2 M N f(Z N ' r *(s))ds 

Jo 

= f(Z N > r *(t)) - f(Z N ^(0)) - [ N r m N f(U 2 Z N ^(s) + y)V^{dy x ds) 

-/Lix[0,t] 

is a martingale. Since / and N ri B>Nf are bounded by constants, N ri ~ T2 M^ is bounded by 
a constant on any bounded time interval. It follows that {N ri ~ r2 M^} is relatively compact, 
any limit point is a martingale with initial value zero, and any limit point is Lipschitz 
continuous with Lipschitz constant sup z |Co/(z)|. Since any continuous martingale with 
finite variation paths is constant, it follows that the limit must be zero. Combining these 
observations with those of the previous section, we have the following theorem. 
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Theorem 5.1 Assume Condition 4-5 Suppose (4-2) holds and that (4-3) is stochastically 
bounded. Let Z T 2 2 and r M be as in the conclusion of Theorem 4-6 Then for all f G C^(R S °), 



C f(Z r 2 2 (s)+y)V 1 (dyxds) 



" Z?(s) f(y)V 1 (dy x ds) = 0. 



LlX[0,Too) 



'LiX[0,roc,) 

If for each z G n 2 E ; tc z is the unique stationary distribution for C z , then 

V x {dy x ds) = Tc z r 2(s) (dy)ds, 
and the system of equations in Theorem \4-6\ becomes 

e ■ z?(t) = e ■ z?(o) + £ (0 ■ ^ Y *( f f x ^ z 2 2 ( s ) + y)* z?{s \d y )ds) 



fcer; 



for 9 G K2 with a$ = 0, and 

6 ■ D a °Z?(t) = 6 ■ D a *Z?(0) + £ (0 ■ ^C*) I I MZ r 2 2 (s) + y)n z ?^(dy)ds, 

, JO Jhl 



fcer; 



for 9 G K2 with a$ > 0, with the equations holding for t G [0,Too)- 

Remark 5.2 Assuming uniqueness, the system determines a piecewise deterministic Markov 
process in the sense of Davis (1993). If one defines 



/ \ k (z + y)ir z (dy), z G n 2 E, 
Jhi 



the description of the system will simplify. 

We still need to address conditions for the relative compactness of the sequence of occu- 
pation measures. If (I — n 2 )E is compact, relative compactness is immediate. Otherwise, it is 
natural to look for some kind of Lyapunov function. Note that if 7^ = inf{t : \Z N,T2 (t)\ > c}, 
then 



f(Z N ^(tA^)-f(Z N ' r2 (0)) 
is a martingale for all locally bounded /. 



tA 7 « 



N 



' - T 



"JV 



f(Z N ^(s))ds 



Lemma 5.3 Let h q fi and ifj q ^ be as in (4-3). Suppose f* e are nonnegative functions and 
there exist positive constants c±, c 2 such that 



sup N r m N f» e (z) <c x - c 2 ip qfi {,h qfi ((l - U 2 )z)) 

N 



for all z satisfying |IT 2 z| < q and for each c£l, 



sup{|(7 - U 2 )z\ : |n 2 z| and sup N r2 M N f" g (z) > c} < 00. 

N 



Then for each t > 0, {V^ ,T2 } is relatively compact and (4-3) is stochastically bounded. 
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6 Examples 



We give some additional examples that demonstrate how identifying exponents satisfying 
the balance condition leads to reasonable approximations to the original model. For a "pro- 
duction level" example, see the analysis of an E-coli heat shock model in Kang (2009). 



6.1 Goutsias's model of regulated transcription 



We consider the following model of transcription regulation introduced in Goutsias (2005) 



and studied further in Macnamara, Burrage, and Sidje (2007). The model involves six species 



and ten reactions: 



X 1 = # of M 
X 2 = #ofD 
X 3 = # of RNA 
X 4 = # of DNA 
X 5 = # of DNA-D 
X 6 = # of DNA-2D 

RNA -> 

M -> 

DNA-D -> 

RNA -> 

DNA + D -> 

DNA-D -> 

DNA-D + D -> 

DNA -2D -> 

M + M -> 



Protein monomer 
Transcription factor 
mRNA 

Unbound DNA 

DNA bound at one site 

DNA bound at two sites 





iZiVA + DiVi 


DNA-D 
DNA + D 
DNA ■ 2D 
DNA -D + D 
D 

2M . 



D 
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Taking the volume V = 1, the corresponding system of equations becomes 

X 1 (t) = X 1 (0) + Y 1 (k[ [ X 3 {s)ds)) + 2Y 10 ( K ' 1Q I X 2 {s)ds) 

Jo Jo 

-Y 2 (k' 2 [ Xi(s)tfe) - 2Y 9 (4 f X^s^X^s) - l)ds 
Jo Jo 

X 2 (t) = X 2 (0) + Y & (k' 6 [ X 5 (s)ds) +r 8 (4 / X 6 (s)ds) + Y 9 (K' g [ - 

</o ./o ./o 

-y 5 («' 5 / X 2 (s)X 4 (s)ds)-Y 7 (K , 7 [ X 2 (s)X 5 (s)ds)-Y w (K' 10 [ X 2 (s)ds) 
Jo Jo Jo 

X 3 (0) + Y 3 (k' 3 [ X 5 (s)ds)-Y i (K , A [ X 3 (s)ds) 
Jo Jo 

X 4 (0) + F 6 k/ X 5 (s)ds)-Y 5 (K' 5 [ X 2 (s)X 4 (s)ds) 
Jo Jo 

X 5 (0)+y 5 («' B / X 2 (s)X 4 (s)ds) + Y 8 (K' s [ X 6 (s)ds) 
Jo Jo 

-y 6 « / X 5 (s)ds)-Y 7 ( K ' 7 [ X 2 (s)X 5 (s)ds) 
Jo Jo 

X 6 (t) = X 6 (0)+y 7 (4 / X 2 (s)X 5 (s)ds)-Y s (4 f Xe(s)ds). 

Jo Jo 

6.2 A scaling with two fast reactions 

In his analysis of the model, Goutsias assumes two time-scales and identifies reactions 9 and 
10 as "fast" reactions. In our approach, that is the same as assuming /? g = /3io > /?i = ■ ■ ■ = 
(3$, so we take N = 100, (3$ = /?io = and 0i = ■ ■ ■ = (3$ = — 1. Recall the relationships 
«4 = k^Nq (we are assuming the volume V = 1) and pk = Pk + v k ■ a. Employing the rate 



*s(f) 
X 4 (t) 

x 5 (t) 



constants from Goutsias (2005), and taking a« = for all i, we have 



Table 1: Scaling exponents for reaction rates 



Rates 


Scaled Rates 


P 




4.30 x 10~ 2 


«i 


4.30 


Pi 


-1 


K 2 


7.00 x 10~ 4 


«2 


0.07 


P2 


-1 


K 3 


7.15 x KT 2 


^3 


7.15 


Pi 


-1 


K,y 


3.90 x 10~ 3 


K 4 


0.390 


P4 


-1 




1.99 x 10~ 2 


K 5 


1.99 


P5 


-1 




4.79 x lO" 1 


K 6 


47.9 


P6 


-1 


K/ij 


1.99 x 10~ 4 


K 7 


0.0199 


P7 


-1 


K 8 


8.77 x 10~ 12 


K 8 


8.77 x 10- 10 


P8 


-1 


K 9 


8.30 x 10~ 2 


Kg 


0.0830 


P9 





K 10 


5.00 x 10- 1 


«10 


0.500 


PlO 
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Then, for 7 = 0, (Z±' , Z 2 ' Q ) converges to the solution of 

Z\{t) = X 1 (0) + 2F 10 Ko / Z°(s)ds)-2Y 9 (K 9 [ Z°(s)(Z°(s) - l)ds 

Jo Jo 

Z° 2 (t) = X 2 (0) + Y 9 (k 9 [ Z°(s){Z°(s)-l)ds-Y 10 ( Kl0 f Z° 2 (s)ds), 

Jo Jo 



and for k > 2, Z^'° converges to X k (0). 



For 7 = 1, the kind of argument employed in (3.12) implies 

-t 

rNA/ \/ryNAt \ ,\J / rzNA 



k 9 / Z^\s){Z^\s) - l)ds - / K 10 Z?'\s)ds ->• (6.1) 
Jo Jo 

but does not lead to a closed system for the limit of {Z^' 1 , . . . , Z^' 1 ). To obtain a closed 
limiting system, we introduce the following auxiliary variable: 

Z%\t) = Z?> 1 (t) + 2Z?> 1 (t) 

= Z^ 2 ' 1 (0)+Y 1 {k 1 I Z?'\s)d S )) + 2Y 6 (K 6 [ Z^ 1 (s)ds) + 2Y 8 (n 8 [ Z^\s)ds) 
Jo Jo Jo 

-2Y 5 (k 5 f Z 2 N \s)Z^\s)ds)-2Y 7 {K 1 f Z^\s)Z^\s)ds)-Y 2 {K 2 f Z^\s)ds) 
Jo Jo Jo 

and observe that the conditional equilibrium distribution satisfies 

Kg(Zi + 2)(Z! + l)fJL 8 (zx + 2,Z 2 -1) + K 10 (z 2 + l)(JL a {Zi - 2, Z 2 + 1) 

= (k S Zi(z! - 1) + K 10 Z2)lI s (zi,Z 2 ) 

and is uniquely determined by the requirement that 

Zi + 2z 2 = Zl 2 (s), 

where Z\ 2 is the limit of Z^ 1 . For m = z\ + 2z 2 , the conditional equilibrium distribution is 

fi m ( Zl , z 2 ) = M m K 107 , 9; , , (6.2) 

Z\\Z 2 \ 

where M m is a normalizing constant making \x m a probability distribution on the collection 
of (zi, z 2 ) such that z\ and z 2 are nonnegative integers satisfying z\ + 2z 2 = m. Define 



mm 



)= I z 2f i m (dz 1 ,dz 2 ) = M m Yl J-° 2 ^L-iV ' (6 ' 3) 

l<z 2 <m/2 V 2J ' ( 2 >' 
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and observe that m — 2a(m) = f Zifi m (dzx, dzi). Then (Z 12 , Z 3 ' , . . . , Z 6 ' ) converges to 
the solution of 

Z l l2 (t) = Z 1 12 (0)+Y 1 (k 1 [ Zl(s)ds)) + 2Y 6 (n 6 [ Zl{s)ds) + 2Y 8 (k 8 [ Zl(s)ds) 

Jo Jo Jo 

-2Y 5 {k 5 [ a(Z l 12 (s))Zl(s)ds)-2Y 7 (K 7 [ a(Z{ 2 (s))Zl(s)ds) 
Jo Jo 

-Y 2 (k 2 I {Z\ 2 {s)-2a{Z\ 2 {s)))ds) 
Jo 

Z\{t) = Zl(p)+Y 3 (K 3 f Z 1 5 (s)ds)-Y 4 (^ [ Zl(s)ds) 

Jo Jo 

Z\(t) = Zl{0)+Y 6 (Kef Z 1 5 (s)ds)-Y 5 { K5 [ a{Z\ 2 {s))Zl{s)ds) 

Jo Jo 

Z\{t) = Z1(0)+Y b (k b [ a(Zl 2 (s))Zl( S )ds)+Y 8 (K 8 [ Zl{s)ds) 

Jo Jo 

-Y 6 (Ke f Zl(s)ds)-Y 7 (K 7 I a{Z\ 2 {s))Zl{s)ds) 
Jo Jo 

Z\{t) = Z^0)+Y 7 {Krf a (Z 1 12 (s))Zl(s)ds)~Y 8 (K 8 [ Zfo)ds) , 

Jo Jo 

which is essentially the approximation obtained by Goutsias. Note that the "fast" reactions, 
reactions 9 and 10, have been eliminated from the model. 

This system is not entirely satisfactory as a(m) is not computable analytically. For 



simulations, values of a{m) could be precomputed using (6.3). E, Liu, and Vanden-Eijnden 



(2007) suggest a Monte Carlo approach for computing a{m) as needed. Goutsias suggests a 



way of approximating the transition rates which is equivalent to the following: The limit in 



(6.1 ) implies 

«i a(m) = K g I zi(zi - l)fi m (dzi,dz 2 ) (6.4) 



as can be verified directly from the definition of fi m . A moment closure argument suggests 



replacing (6.4) by 

K W a(m) = k 9 I z 1 /j, m (dzi,dz 2 ) I (z 1 - l)/j, m (dzi, dz 2 ) 



= Kg(m — 2a(m))(m — 2a(m) — 1), 
which gives a quadratic equation for the approximation for a(m). 

6.3 Alternative scaling 

Observe that k' 9 < k' 6 , so reaction 6 is actually "faster" than reaction 9. Consequently, it 
is reasonable to look for a different solution of the balance conditions with /?io = > f3g. 
Drop the assumption that «j = 0, and consider a subset of the balance equations. Recall 
that p k = /3 k + v k ■ a. 
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Table 2: Balance equations 



Variable 


Balance equation 


x 1 


Pi V P10 = P2 V p 9 


x 2 


P6 V p 8 V p 9 = p 5 V p 7 V pio 


x 3 


P3 = P4 


x 4 


P5 = P6 


x 5 


P5 V p 8 = p 6 V p 7 


x & 


P7 = P8 


X l + 2X 2 + 2X 5 + 4X 6 


Pi = P2 


X 2 + X 5 + 2X 6 


P9 = PlO 


^5 + Xq 


P5 = P6 


X4 + X 5 + Xq 


= 


X4 + x 5 


P8 = P7 



We take No = 100, cti = a 2 = 1, and = for 3 < i < 6. We see that the following 
exponents satisfy the balance conditions and the additional requirement that K k > k\ imply 
fik > A) except for /3 8 , the exponent associated with the extremely small rate constant k' 8 . 
Recall that n k is determined by the requirement n' k = K k N^ k . 



Table 3: Scaling exponents for reaction rates 



Rates 


Exponents 


Scaled Rates 


P 


K\ 


4.30 x 10~ 2 


Pi 


-1 




4.30 


Pi 


-1 


K 2 


7.00 x 10~ 4 




-2 


«2 


7.00 


P2 


-1 


K 3 


7.15 x 10~ 2 




-1 


«3 


7.15 


P3 


-1 




3.90 x 10- 3 


Pa 


-1 


K 4 


0.390 


Pa 


-1 




1.99 x 10~ 2 


Pb 


-1 


K 5 


1.99 


P5 





Kg 


4.79 x 10" 1 


Pe 





Kq 


.479 


P6 







1.99 x 10~ 4 


^7 


-3 


K 7 


199 


P7 


-2 


Kg 


8.77 x 10~ 12 


Ps 


-2 


Kg 


8.77 x 10~ 8 


P8 


-2 


Kg 


8.30 x 10~ 2 


^9 


-1 


Kg 


8.30 


P9 


1 


K 10 


5.00 x 10- 1 


PlO 





«10 


0.500 


PlO 


1 
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Defining Z? ,7 (t) = N^X^(Nn) and n k = N~^K' k , 

Z?'\t) = Z? ,7 (0) + N-'Y^ f ^N^Z^is) ds) + 2iV- 1 r 1 o( f k 10 W +1 Z?" (s) ds) 

Jo Jo 

-N- l Y 2 { [ ^N^Z^'is) ds) - 2N- 1 Y 9 ( [ k 9 N^ +1 Z? • 7 (s)(Z 1 Ar ' 7 (s) - N' 1 ) ds) 
Jo Jo 

Z?«(t) = Z^(0) + iV- 1 F 6 ( f KeWZ^is) ds) + N'^ f k 8 N^ 2 Z^(s) ds) 

Jo Jo 

+N' 1 Y 9 ( f k 9 N^ +1 Z^(s)(Z^(s) - N" 1 ) ds) - N-'Ysi f k 5 N~<Z^(s)Z^(s) ds) 
Jo Jo 

-N^Yri f k 7 N^- 2 Z^(s)Z^(s) ds) - ^ 1 F 10 ( l\ w N^ +1 Z^{s) ds) 
Jo Jo 

Z^(t) = Z^(0) +Y 3 ([ t KsN^Z^is) ds) - Y 4 ( f k,N^ 1 Z^(s) ds) 

Jo Jo 

Zf\t) = Zf ' 7 (0) + y 6 ( f KeWZ^(s) ds) -Y 5 {f k 5 WZ^(s)Z^(s) ds) 

Jo Jo 

Z^(t) = Zf' 7 (0) + Y 5 ( f k 5 N~<Z?>\s)Z^(s) ds) + Y 8 ( f k 8 N^ 2 Z^\s) ds) 

Jo Jo 

-Y 6 ( f KeWZ?"( 8 ) ds) - Y 7 ( f k 7 N^ 2 Z^\s)Z^'{s) ds) 
Jo Jo 

Zf"(t) = Z^(0) + Y 7 ( f k 7 W- 2 Z^{s)Z^{s) ds) - Y 8 ( f k 8 N^ 2 Z^(s) ds). 

Jo Jo 

Useful auxiliary variables include 

NZ? ,7 (t) + 2NZ^(t) + 2Z^(t) + 4Z^(t) = NZ?(0) + 2NZ 2 N (0) + 2Zf(0) + 4Z? (0) 

+Y 1 ( [ ^N^Z^is) ds) - Y 2 ( [ k 2 N^ 1 Z*'\s) ds) 
Jo Jo 

NZ^(t) + Z^{t) + 2Z^(t) = NZ?(0) + (0) + 2Z* r (0) 

+Y 9 ( [ k 9 N^ +1 Z^(s)(Z^(s) - N- 1 ) ds) - Y 10 ( [ k w N^ +1 Z^{s) ds) 
Jo Jo 

Z»"(t) + Z^(t) = Z»(0) + Z^(0) + Y 5 ( J' k 5 N^Z^(s)Z^(s) ds) - F 6 ( J* k 6 N^Z^(s) ds) 
Z»"(t) + Z^(t) + Z^(t) = Z» (0) + Z» (0) + Z»(0) 

Z»"(t) + Z^(t) = Z» (0) + Z» (0) + Y 8 ( f k 8 N~<- 2 Z%<\s) ds) - Y 7 ( f k 7 N^ 2 Z^(s)Z^(s) ds) 

Jo Jo 
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For 7 = 0, the limiting system is the piecewise deterministic model 



Z°(t) 


= Z 1 (0)+ [ 
Jo 




» - 2/? 9 Z 1 °(s) 2 ) ds 




z° 2 (t) 


= Z 2 (0) + / 
Jo 


Mis] 


) 2 - K 10 Z°( S )) ds 


(6.5) 


z!(t) 


= Z 4 (0) + F 6 ( 


f « 6 Z 5 °l 
Jo 


[s)ds)-Y B (f k 5 Z° 2 (s)Z° 4 (s) ds) 
Jo 




zm 


= Z 5 (0) + Y 5 ( 


[ k 5 Z 2 °i 
Jo 


[s)Z° A (s)ds)-Y 6 ([ K 6 Z° 5 (s)ds), 
Jo 





with Z$(t) = Z 3 (0) and Z°(t) = Z 6 (0). 

For 7 = 1, we introduce the auxiliary variables 



Z N /(t) = ^ 1 (t) + 2^ 1 (t) 



= Zf(0) + Zf(0) + r 8 (/ K 8 N- 1 Z^ 1 (s)ds)-Y 7 ([ K 7 N- 1 Z 2 N '\s)Z^ 1 (s)ds). 

Jo Jo 

Observing that Z^' 1 is asymptotically the same as Z^' 1 + 2Z 2 ' 1 + 2N~ 1 Z^ 1 + AN _1 Zq' 1 , 
Z^' 1 converges to Z\ 2 (t) = Zi 2 (0) = limjv->oo(Zf r (0) + 2Z^ v (0)). In particular, Zj 2 is constant 
in time. We also have Z 4 1 5 (t) = Z 45 (0) = lim^ oo (Z 4 Ar (0) + Z^(0)). 

Let l/^' 1 denote the occupation measure for {Z x ' , Z 2 ' , Z 4 ' , Z 5 ' ). The stochastic 
boundedness of Z^' 1 and Z^' 1 ensures the relative compactness of {V^' 1 }, and as in Section 
5| V 1 ^' 1 converges to V 1 (dz,ds) = v s (dz)ds where v s satisfies 

Cfv s (dz) = 
and 

Cf(zi,z 2 ,z 4 ,z 5 ) = (k w z 2 - K 9 zf)(2d z J(z) - d Z2 f{z)) 

+K 6 z 5 (f(z + e 4 - e 5 ) - f(z)) 
+n b z 2 z A (f(z - e 4 + e 5 ) - f{z)). 

Consequently, v s is uniquely determined for each s by the requirement that z\ + 2z 2 = 
Z\ 2 {s) = Z 12 (0) and z 4 + z 5 = Z 45 (s) = Z 45 (0), and hence, 



v s (dz) = 5 vli z 12 (o))(dz 1 )5 V2{Zl2m (dz 2 )B(Z 45 (0), /n w > dz ^ dz 5) 



^6 

«6 + ^5V ? 2(Zi 2 (0))' 



where 

a/ Kf + 8K 9 K 10 y - KlQ 



My) 
My) 



4k 9 

4k 9 ?/ + K 10 - a/ k 2 10 + 8K 9 Kioy 
8 Kg 
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and B( n } p, dz 4} dz§ ) is given by the binomial distribution 



P{Zt = k,Z B = n-k}= ( 



Averaging gives 



(6.6) 



Zj(t) = Z 8 (0) + F 3 ( T ^5^(^(0)) 45(Q) ds) _ y4( r ^ 1(s) dsy 

N 2 

Finally, for 7 = 2, dividing the equation for Z 3 ' by iV, we see that 

f Q ZF{s)dsn^f Q Z»*(s)d8 t 

and (Z 12 2 , Z^ 2 , Zq' 2 ) converges to the solution of 

ZUt) = Z 12 (0) + J* (^Zl(s) - wi(ZUs))) ds 

ZUt) = Z A5 (0) + Y s ([ K8 Z 2 (s)ds)-Y 7 ([ K7V2 (Z 2 2 (s))zl(s)ds) 

Jo Jo 

Zlit) = Z 6 (0) + Y 7 ([ K 7 <f 2 (Z 2 12 (s))zl(s)ds)-Y 8 ([ K 8 Z 2 (s)ds) 

Jo Jo 

7 2 (+\ - K ^2(Zl 2 (t)) v2 

K 6 + K 5 ip 2 (Z( 2 {t)) 

6.3.1 Simulation results 

We compare simulation results for the full model with the approximations given by the 
limiting systems. The mean and standard deviations of the number of molecules for each 
species or for the auxiliary variables of interest are given from 100 simulations of the full 
model and from 1000 simulations of the limiting systems. The evolution of the processes in 
the full model is approximated by the evolution of the processes in the limiting system using 
the relationship 

Xi(t) = X^{t) « N^Z]{tNp). 



Following Goutsias (2005), initial values are taken as Xl(0) = 2, X 2 (0) = 6, X 5 (0) = 2, and 
all other values equal to zero. 

For 7 = 0, we observe the evolution of the processes during the time interval [0, 100]. 



The full model is reduced to the 4-dimensional hybrid model (6.5) in which Z\ and Z® are 
the solution of a pair of ordinary differential equations and Z\ and Z® are discrete with 
transition intensities depending on Z\. The evolution of X\, X 2 , X 4 , and X 5 in the full 
model is given in Figure [T] and the evolution of the approximation is given in Figure [2] The 
exact simulations of the full model are done using Gillespie's stochastic simulation algorithm 



(SSA) from Gillespie (1977). For the approximation, Z\ and Z® are solved by the Matlab 
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Figure 1: Simulation of the full model during t = to t = 100 
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Figure 2: Approximation using the limiting model for 7 = in the alternative scaling 
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ODE solver, and Z\ and Z\ are computed by Gillespie's SSA taking Z 2 from the solution 
of ODE. The evolution of X\ and X 2 are well captured by Z\ and Z 2 in Figure [2j These 
deterministic values approximate the evolution of the mean of X\ and X 2 given in Figure [T] 
except for a slight increase over time in the simulation of the full model. Note that in the 
approximate model Z®(t) + 2Z 2 (t) is constant, but that is not the case in the full model. 
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Figure 3: Simulation of the full model during t — to t — 1000 
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Figure 4: Approximation using the limiting model for 7 = 1 in the alternative scaling 




For 7 = 1, we consider the evolution of the processes on the time interval [0, 1000]. The 



full model is reduced to the 1-dimensional limiting system (6.6 ) with a single jump process Z\ . 
Comparing the governing equations for Z^' 1 and Z\ , the different behavior of the evolution of 

the two processes comes from the difference between Z^' 1 and Z 5 (t) = ^+^2(2^2(0)) ^ 45 W • 
Therefore, plots of the evolution of both X 3 and X 5 in the exact simulation are given in 



Figure 3 In Figure [4J the evolution of Z\ and o 
simulations, we use Gillespie's SSA. In Figure 



Z 5 is given. For both exact and approximate 
3j Z5' 1 increases slightly and then decreases 
the increase during the early time 



to zero. Since Z 5 is approximated as constant in Figure 
and the decrease to zero of X 3 is not captured by the approximation. 

For 7 = 2, the simulation is carried out on the time interval [0, 10000]. The 3-dimensional 
limiting model (6.7) is piecewise deterministic and includes the auxiliary variables Z\ 2 



7 2 

^45' 



and the species abundance Z\. Z\ 2 is governed by a random differential equation driven by 
a component of the jump process, Z| 5 . Z\ h and Z| are discrete with transition intensities 
that depend on Z\ 2 . Since there is mutual dependence between the continuous and discrete 
components, we modify Gillespie's SSA to simulate the limiting system. Here is a brief 
description of the simulation method for the limiting system. 



1. Assume that the process has been simulated up to ti, the zth. jump time of the jump 
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Figure 5: Simulation of the full model during t = to t = 10000 
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Figure 6: Approximation using the limiting model for 7 = 2 in the alternative scaling 
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process. Simulate a unit exponential random variable A by simulating a uniform [0, 1] 
random number r 1 and setting A = log 

2. Solve the differential equation for Z\ 2 starting at Zf 2 (ti) holding Z% 5 (t) = Z\ h {ti) and 
Z\{t) = Zg(ti) until time t i+ \ satisfying 



rti+i _ 2 
/ (K 7 cp 2 {Zl 2 (s))Z 5 (s) + K 8 Z%(s))ds 

r u ^ K 5 K 7 <p 2 (zUs) 

l t . K 6 + K 5 <f 2 (Zf 2 (s)) 



Z%^i) I 5 7^2 ( 12 _ i _|_ Ks Z£(ti)(t i+ i - ti 

Ju K 6 + K 5 (fi2(Zi 2 (s)) 



= A. 

(We compute the integral by the trapezoid rule using the grid points from the ODE 
solver.) 

3. Simulate a uniform [0, 1] random number r 2 . If 

2 

r2 < K 7 <p 2 (Zf 2 (ti + i))Z 5 (t i+1 -) , g7 , 

K 7 (f 2 (Zf 2 (t i+1 ))Z 5 (t i+1 -) + K 8 Zl{t i+l -)) 

K 5 K 7 y 2 (Zf 2 (t i+1 )) 2 Zl 5 (ti) 

K 5 K 7 Cp 2 (Zf 2 (ti +1 )) 2 Zl 5 (ti) + « 8 Z|(tj)(« 6 + K h (f 2 {Zl 2 {t i+1 ))y 

set 

zUu +l ) \ _ ( zuu) \ , ( -i 



and if the reverse inequality holds in (6.7), set 



4. Go back to step 1. 

Comparing plots for Xi(t) + 2X 2 (t) in Figure [5] and for N Z 2 2 (tNQ 2 ) in Figure |6j the plot 
in the approximation increases more rapidly at early times and starts to drop earlier than 
the plot in the exact simulation. Also, the peak level in the approximation is much lower 
than the peak level in the exact simulation. 

Since k 8 = 8.77 x 10~ 8 is small compared to the time interval, Reaction 8 will rarely 
occur on the time scales we are considering. We retained this reaction in the limiting model 
only to emphasize that a long time after the model appears to equilibrate, action may restart 
after the dissociation 

DNA ■ 2D ->> DNA ■ D + D. 

If Reaction 8 does not occur, the stochastic behavior of the limiting model just depends 
on the two jump times 

r 2 = inf{t : Z^it) = 1}, r 2 = M{t : Z 2 45 (t) = 0}, 
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so we compare these random variable to the corresponding variables 



n = inf {t : X 4 (t) + X B (t) = 1}, r = inf {t : X 4 (t) + X 5 (t) = 0}, 

from the original model or more precisely, because of the change of time scale, we compare 
(iVgr?,jy?7f)to(r X) 75). 

In Figure |5j plots for t\ and To for 100 exact simulations are given. Taking the average, 
the mean of first hitting time of X 4 (t) + X 5 (t) to 1 is 305.44 and the mean of the first hitting 
time of X 4 (t) + X 5 (t) to is 512.45. In Figure |6j plots for 1000 simulations of rf and Tq 
are given. The mean of the first hitting time of Z 2 5 (tN Q ~ 2 ) to 1 is 155.95 and the mean of 
the first hitting time of Z^tN^ 2 ) to is 261.01. Comparing the two stopping times in the 
simulations of the full model and of the approximation, the mean hitting time to 1 and in 
the approximation is much faster than in the full model. Consequently, the quicker decrease 
of Z| 5 to gives a discrepancy in the peak levels and the peak times in the full model and 
in the approximation. 

6.4 Derivation of Michaelis-Menten equation 



Darden (1979, 1982) derives the Michaelis-Menten equation from a stochastic reaction net- 
work model. His result can be obtained as a special case of the methods developed here. 
Consider the reaction system 

K 'i k' 

S% + S 2 ^ S3 ^ S 4 + S 2 , 

K 2 

where Si is the substrate, S 2 the enzyme, S3 the enzyme-substrate complex, and S 4 the 
product. Assume that the parameters scale so that 

= Z^(Q)-N- l Y x {N f 1 K 1 Z^(s)Z^(s)ds) + N- l Y 2 {N f \ 2 Z${s)ds) 

Jo Jo 

Z»{t) = Zf(0)-n(iv/ Kl Z»{s)Z?{s)ds)+Y 2 {N [ K 2 Z£{s)ds) 

Jo Jo 

+Y 3 (N [ K 3 Z?(s)ds) 
Jo 

Z?(t) = Z 2 N (0) + Y 1 (N [ K 1 Z^(s)Z 2 N (s)ds)-Y 2 (N [ k 2 Z* \s)ds) 

Jo Jo 

-Y 3 (N [ K3 Z?(s)ds) 
Jo 

Z»{t) = N-^N [ n 3 Z?(s)ds), 

Jo 

that is, a\ = a 4 = 1, a 2 = a 3 = 0, (3\ = 0, and (3 2 = /3 3 = 1. 
Note that M = Z 3 (t) + Z 2 {t) is constant, and define 



V 2 N {t)= f Z»{s)ds. 
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Theorem 6.1 Assume that Z^(0) — > £a(0). Then (Z 1 ^ ,V 2 N ) converges to (x\(t),v 2 (t)) 
satisfying 

xi(t) = xi(0)- [ K 1 xi{s)v 2 {s)ds + [ k 2 (M - v 2 (s))ds (6.8) 
Jo Jo 

— — Kix 1 (s)v 2 (s)ds + / (k 2 + k 3 )(M — v 2 (s))ds, 
Jo Jo 



and hence v 2 (s) = — , , — ^ and 



(k 2 +K 3 )M 

MniK 3 Xi(t) 



Xl(t) 



K 2 + K 3 + KiXAs) 



Proof. Relative compactness of the sequence (Z± , V 2 ) is straightforward. Dividing the 
second equation by N and passing to the limit, we see that any limit point (xx,v 2 ) of 
(Zf, V 2 N ) must satisfy 

= - / K 1 xi(s)dv 2 (s) + (k 2 + n 3 )Mt - / (« 2 + K 3 )dv 2 (s). (6.9) 
Jo Jo 



Since v 2 is Lipschitz, it is absolutely continuous, and rewriting (6.9) in terms of the derivative 



gives the second equation in (6.8). The first equation follows by a similar argument. □ 



6.5 Limiting models when the balance conditions fail 



The balance condition, Condition 3.2 has as its goal ensuring that the normalized species 



numbers remain positive, at least on average, and bounded. Mastny, Haseltine, and Rawlings 



(2007) consider examples in which model reduction is achieved by eliminating species whose 
numbers are zero most of the time. We translate some of their examples into our notation 
and see how one can obtain reduced models even though the balance conditions fail. 
Consider 

K ' K , 
Si^j2S 2 , S 2 ^S 3 , 

where we assume k 2 , k' 3 » k[. We take the scaled system to be 

Z?(t) = Z 1 (0)-Y 1 ([ Kl Z?(s)ds)+Y 2 (N [ K 2 Z^{s){Z^{s)-l)ds) 

Jo Jo 

Z»(t) = Z 2 (0) + 2Y 1 ([ Kl Z»{s)ds) - 2Y 2 (N [ k 2 Z* '(s)(Z? '(s) -l)ds) -Y 3 {N [ n 3 Z 2 N (s)ds) 

Jo Jo Jo 

Z»(t) = Z 3 (0) + Y 3 (N [ K 3 Z 2 N {s)ds). 

Jo 

Consequently, assuming Z?(Q) = 0, for most t > 0, Z?(t) = and 

2Y 1 ([ K 1 Z?(s)ds)=Y 3 (N [ K 3 Z 2 N (s)ds) + 2Y 2 (N [ k 2 Z* \s){Z* \s) - l)ds) . 
Jo Jo Jo 
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To be precise, letting A denote Lebesgue measure and defining 



for each t > 0, 



l {Z»(r-)=2] 



dR»(r), R»(t) 



Hm A{0<s<t:Z?(s)^0}< lim / Z£ (s)ds 

7V-»oc N-^-oo 







lim sup sup Z 2 (s) < 2 

N^oo s<t 



lim / \R" (a) - R% (s)\ds = 



lim / \R$(s) -2R 3 y (s)\ds = 



so 



lim / \R? (s) - R% (s) - R$ (s)\ds = 0. 



Setting Q N (t) = 1 



R%{t)- I NQ N {s)K 2 2ds and R$ (t) - I NQ N {s)K 3 2ds 



are martingales. Working first with a subsequence satisfying (A. 7), by Lemma A. 13, (R 2 , R 3 ) 
converges to counting processes (R2, R3) with intensities 



K2 + K 3 



«2 + «3 



where Z\{t) = Z\(0) — R 3 (t). It follows that the finite dimensional distributions of (Z± , Z 3 ) 



converge to those of a solution to 

Z x {t) = Z 1 (0)-Y([ t 

Jo 



Z 3 (t) = Z 3 (0) + 2Y(f t 

Jo 



^1^3 

«2 + ^3 
«4«3 



^2 + ^3 



Zi(s)ds) 
Z 1 (s)ds) 



which is the reduced model obtained in Mastny et al. (2007). More precisely, (Z^,Z. N " 



A.14 



converges in the Jakubowski topology as described in Remark 

(Note the relationship between our rate constants and those of Mastny et al. (2007): 
K\ = fci, k 2 = \k-x, and k 3 = k 2 .) 
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A Appendix 

A.l Convergence of random measures 



The material in this section is taken from Kurtz (1992). Proofs of the results can be found 
there. 

Let (L, d) be a complete, separable metric space, and let Ai(h) be the space of finite 
measures on L with the weak topology. The Prohorov metric on .M(L) is defined by 

p(/i, v) = inf {e > : p(B) < u(B 6 ) + e, u(B) < p(B e ) + e, B G B(L)}, (A.l) 

where B e = {x G L : inf yE s d(x, y) < e}. The following lemma is a simple consequence of 
Prohorov's theorem. 

Lemma A.l Let {r n } be a sequence of M.{h) -valued random variables. Then T n is rela- 
tively compact if and only if {r n (L)} is relatively compact as a family of M.-valued random 
variables and for each e > 0, there exists a compact K C L such that sup n P{r n (f^ c ) > e} < 
e. 

Corollary A.2 Let {r n } be a sequence of j\4(h) -valued random variables. Suppose that 
sup n -E[r n (L)] < oo and that for each e > 0, there exists a compact K C L such that 

limsup J E[r n (K c )] < e. 

Then {T n } is relatively compact. 

Let £(L) be the space of measures on L x [0, oo) such that /x(L x [0,t]) < oo for each 
t > 0, and let £ m (L) C £(L) be the subspace on which /x(L x [0, £]) = t. For fi G £(L), let /x* 
denote the restriction of p to L x [0, t). Let pt denote the Prohorov metric on M(L x [0, t}), 
and define p on £(L) by 



oo 



that is, {/i n } converges in p if and only if {p^} converges weakly for almost every t. In 
particular, if p(p n , p) — >■ 0, then Pt(Pn, p l ) — > if and only if p„(L x [0, t]) —> p{L x [0, t]). 



The following lemma is an immediate consequence of Lemma A.l 



Lemma A. 3 A sequence of (£ m (L),p) -valued random variables {T n } is relatively compact 
if and only if for each e > and each t > 0, there exists a compact K C L such that 
mi n E[T n {Kx [0,t])} > (l-e)t. 

Lemma A. 4 Let T be an (£(L),p) -valued random variable adapted to a complete filtration 
{Ft} i n the sense that for each t > and H G B{h), T(H x [0,t]) is j r t -measurable. Let 
\(G) = T(L x G). Then there exists an {Ft} -optional, V{L)-valued process 7 such that 



h(y,s)T(dyxds)= h(y,s) ls (dy)X(ds) (A.2) 

Lx[0,t] JO </L 

for all h G B(h x [0, 00)) with probability one. If A([0,t]) is continuous, then 7 can be taken 
to be {F} -predictable. 
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Lemma A. 5 Let {(x n ,/i n )} C D K [0, oo) x £(L), and (x n , fi n ) — > (x,fi). Let h G C(E x L) 
and if) be a nonnegative function on [0, oo) satisfying lim^oo ip(r)/r = oo sitc/i £/ia£ 



sup / ip(\h(x n (s),y)\)fi n (dy x ds) < oo (A. 3) 

n Jhx [0,t] 

for each t > 0. 

u n (t) = / h(x n (s),y)n n (dy x ds), u(t) = / h(x(s),y)n(dy x ds) 

Jhx[0,t] Jhx[0,t] 

z n (t) = fi n (h x [0,t]), and z(t) = fi(h x [0,t]). 

a) If x is continuous on [0,t] and lim^oc z n (t) = z(t), then lim^oo u n (t) = u(t). 

b) If (x n , z n , /i n ) -> (x,z,fi) in £> ExM [0, oo) x £(L), then (x n , z n ,u n , /i n ) (x,z,u,fi) in 
-DexMxr[0, oo) x £(L). In particular, lim^oo u n (t) = u(t) at all points of continuity of 
z. 

c) The continuity assumption on h can be replaced by the assumption that h is continuous 
a.e. v t for each t, where v t 6 M(E x L) is the measure determined by v t (A x B) = 
fi{(y,s) : x{s) G A,s <t,y G B}. 



The Lemma A.5 and the continuous mapping theorem give the following. 



Lemma A.6 Sup pose (Z N ,V N ) (Z,V) in D E [0,oo) x £ m (L). Let h G C(E x L) and 



■ip be as in Lemma A.5 If {J * ip(\h(Z N \s),y)\)V N (dy x ds)}is stochastically bounded for all 



t > 0, then 

h(Z N (s),y)V N (dyxds)^ [ h(Z (s) , y))V (dy x ds). 

Lx[0,-] Jhx[0,-] 

A. 2 Martingale properties of counting processes 

A cadlag stochastic process R is a counting process if R(0) = and R is constant except for 
jumps of plus one. If R is adapted to a filtration {Ft}, then a nonnegative {J-" t }-adapted 
process A is an {J^}- intensity for R if 

M{t) = R{t) - [ \{s)ds 
Jo 

is an {J^j-local martingale. Specifically, letting n denote the Zth jump time of R, 

M T '(t) = M(t A n) = R(t An) - I \{s)ds 

is an {J-^j-martingale for each I. 

For simplicity, we assume that A is cadlag. 



o 
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Remark A. 7 For Rf. defined in (2.1) and {Ft} = (t(Ri(s) : s < t,l — l,...,ro), the 
intensity for R k is t — >■ Afc(X(£)). 

Lemma A. 8 For each t > and eaca Z ; 

/•tATj 

I > E[R(t A n)] = E[ X(s)ds] (A.4) 

and 

E[R(t)) =E[f \{s)ds), 
Jo 

where we allow oo = oo. If E[R(t)] < oo for all t > 0, t/ien 

i?(t)- / A(s)ds 
Jo 

is an {Ft} -martingale. 

Two counting processes, R%, R2, are orthogonal if they have no simultaneous jumps. 

Lemma A. 9 Let Ri, . . . ,R m be pairwise orthogonal {Ft}- adapted counting processes with 
{Ft} -intensities X k . Then, perhaps on a larger probability space, there exist independent unit 
Poisson processes Y\, . . . , Y m such that 



R k (t) = Y k ([ X k (s)ds), 
Jo 



and R = J2T=i Rk ^ s a counting process with intensity A = J2T=i ^fc- 
If T\ is the Ith jump time of R, then 



P{R k (n) - R k (n-) = l\F n } = (A.5) 



Remark A. 10 Note that the right side of ( A.5) involves the left limits of the intensities. If 
the intensities are not cadlag, then Afc(rj— ) should be replaced by 

limsup/i -1 / Xk(s)ds. 

ft->0+ Jri-h 

The intensity of a counting process does not necessarily uniquely determined its distribu- 
tion. For example, consider the system 

Ri{t) = Y x {l X{R x {s))ds) 



Rz(t) = Y 2 ([ A(i2i(s))ds). 
Jo 



The intensity for each component is X(Ri(t)), but the two components will not have the same 
distribution. 
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Proof. See iMeyerl (19711) and iKurtzl fll980b . 



□ 



Lemma A. 11 Suppose that R±, . . . , R^ are pairwise orthogonal counting processes adapted 
to a filtration {J 7 ^} with {J 7 ^} -intensities Xf , . . . , A^. Let Aj? (t) = f Q X% (s)ds, and sup- 
pose that (A^,...,A^) =>■ (A 1; ...,A m ) in the Skorohod topology on Z} R m[0, oo). Then 
{(Ri , • • • , Rm)} i s relatively compact in the Skorohod topology and any limit point (Ri, ■ ■ ■ , R m 
consists of pairwise orthogonal counting processes. 
At least along a further subsequence, 



(Aj , . . . , A m , R l , . . . , R m ) =^ (Ai, . . . , A m , Ri, ... , R m ), 

and letting {J 7 ^^} be the filtration generated by (Ai,...,A m ,Ri,...,R m ), Rk — A k are 
{J 7 ^' R }-local martingales and there exist independent unit Poisson processes (Y[, . . . ,Y m ) 
such that 

R k (t) = Y k (A k (t)), k = l,...,m. (A.6) 



Remark A. 12 If the A k are adapted to {J 7 ^}, then R will be the unique solution of (A.6) 
and R R in the Skorohod topology. 



Proof. See Kabanov, Liptser, and Shiryaev (1984) 



□ 



In Section 6.5, we consider an example for which the integrated intensities did not have 
a continuous limit. The next lemma covers that situation. 



Lemma A. 13 Suppose that Rq , R± , . . . , R^ are counting processes adapted to a filtration 
{J 7 ^}, and Ri , . . . ,R% are pairwise orthogonal. Suppose Rq has {J 7 ^} -intensity Xq , and 
Ri, . . . , R% have {J 7 ^} -intensities = NQ N ^ , where Q N > 0. Suppose 



{ \ N N 1\ \ ( \ \ 

[A Q , Hi , . . . , fi m ) =^ {*0, Ml J • • • ) f^m), 



.N 



and 



k=l 



(A.7) 



(A.8) 



for each t > 0. Then {(Rq,R^,. . . , R^)} is relatively compact in the Jakubowski topology 
and for any limit point (Rq, Ri, . . . , R m ), 



R(\ 



k=l 



and Ri, . . . , R m are pairwise orthogonal counting processes with intensities 

Hk(t) 



X k (t) 



A (t). 
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Remark A. 14 The sequence may not be relatively compact in the Skorohod topology since 
we have not ruled out the possibility that the sequence has discontinuities that coalesce. See 



the example in Section 6.5 



The Meyer-Zheng conditions (Meyer and Zheng (1984)) imply relative compactness in 
the Jakubowski topology (Jakubowski (1991)). A sequence of cadlag functions {x n } converges 



to a cadlag function x in the Jakubowski topology if and only if there exists a sequence of 
time changes {^ n } such that (x n o 7„,7„) — >■ (x o 7,7) in the Skorohod topology. ( See \Kurtz 



(1991).) The time-changes are continuous, nondecreasing mappings from [0, 00) onto [0, 00) 
but are not necessarily strictly increasing. Convergence implies C \x n {s) — x(s)\ A Ids — > 0. 
In contrast to the Skorohod topology, if X fi ' X and y n — > y in the Jakubowski topology, then 
niVn) (x,y) in the Jakubowski topology on cadlag functions in the product space. 



x 



Proof. By Lemma A. 11 {Rq} is relatively compact in the Skorohod topology and hence 
in the Jakubowski topology. Let 

m 

K o - / ; K k ■ 

k=l 



The stochastic boundedness of {Rq (t)} for each t > and (A. 8) imply the stochastic bound- 
edness of {Rq (t)} for each t > which by (A. 4) implies the stochastic boundedness of 



Let 7tv be defined by 



{ NQ N {s)Y,^{s)ds}. 
Jo k=i 

p N (t) rn 

/ (l + NQ N (s)J2^(s))ds = t. 
Jo k=i 



Since |7at(s) — 7jv(0I — I s ~ {7iv} * s relatively compact. Define 



\ k (s)ds, 



and observe that 



Af°7iv(t) 



NQ N o lN (s)fMf f o lN (s) 



-ds 



, l + NQ N o lN ( s )j: k ^o lN (s) 
is also Lipschitz with Lipschitz constant 1. Since {7Ar(t),i > 0} are stopping times, 

irf-Afo 7iV 

are martingales with respect to the filtration {J 7 ^^}. 

The Lipschitz properties imply the relative compactness of {(A^ o 7^, . . . , o 7^, 7^)} 
in the Skorohod topology which in turn, by Lemma A.ll[ implies the relative compactness 
of 

{(Af o 77V , . . . , A^ o 7Ar , 7iV) Ri ° 7jv, • • • , RZ In)}- 
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Relative compactness of this sequence in the Skorohod topology ensures relative compactness 
of {(Ri , ■ ■ ■ ,Rm)} i n the Jakubowski topology which in turn implies relative compactness 
of {(Rq , Ri , • • • , Rm)} i n the Jakubowski topology 

Along an appropriate subsequence, we have convergence of 7 a? to a limit 7, 



' NQ N o lN {s)Y. k ^ N k o lN {s) 



convergence of to 



l + NQ"o lN ( s )Z k K°lN(s 



Jo Ei^i°7(s) 



and convergence of (Rq , -Rf , • ■ • , R%) i n t ne Jakubowski topology to a process satisfying 

TO 

Rq = 

fc=l 

Since _R 7(t) — Jq 7 ^ A (s)c?s is a martingale, we must have 

p(t) 

/ \ {s)ds = A(t) 
Jo 

and 

A*(*) = f Jt k ° l{S ) ^ o 7(sW(s)ds. 



EiW°7(s) 

Since i?o is a counting process, the R k , k = 1, . . . , m, must be orthogonal, and i?fc must have 
intensity j^j^o- D 
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